{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# IEEE European Low Voltage Test Feeder: \n",
    "\n",
    "http://sites.ieee.org/pes-testfeeders/resources/\n",
    "\n",
    "The current IEEE test cases are focused on North American style systems; however it is common outside of North America to see low-voltage distribution systems, both radial and meshed. It is important to make sure that tools support both dominant styles of distribution system configuration. This test case seeks to fill a benchmark gap by presenting a number of common low-voltage configurations. This circuit also introduces quasi-static time series simulations.\n",
    "\n",
    "IEEE European LV network is a generic 0.416 kV network serviced by one 0.8 MVA MV/LV transformer and a 11kV external grid. The network supplies 906 LV buses and 55 1-PH loads.\n",
    "\n",
    "\n",
    "# Time Series data from Benchmark dataset\n",
    "In this example the time series data from benchmark is used to perform a time series load flow analysis\n",
    "in pandapower\n",
    "\n",
    "The network data from csv is used to create Powerfactory and pandapower networks.\n",
    "\n",
    "The Sum total of power and maximum voltage is plotted with respect to 1440 time steps( 1 min interval) during the day:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pandas as pd\n",
    "import tempfile\n",
    "from pandapower.pf.runpp_3ph import runpp_3ph\n",
    "from pandapower.control import ConstControl\n",
    "from pandapower.timeseries import DFData\n",
    "from pandapower.timeseries import OutputWriter\n",
    "from pandapower.timeseries.run_time_series import run_timeseries\n",
    "from pandapower.networks import ieee_european_lv_asymmetric\n",
    "from pandapower import pp_dir\n",
    "import warnings\n",
    "warnings.simplefilter(action='ignore', category=FutureWarning)\n",
    "\n",
    "parent = os.path.join(pp_dir,'test','test_files')\n",
    "path = os.path.join(parent, \"European_LV_CSV\")\n",
    "load_path = os.path.join(path, \"Load_Profiles\")\n",
    "def remove_comments(f):\n",
    "    '''Pass comments'''\n",
    "    start=f.seek(0)\n",
    "    for index in range(5):\n",
    "        start=f.tell()\n",
    "        if f.readline().startswith('#'):\n",
    "            continue\n",
    "        else:\n",
    "            break      \n",
    "    f.seek(start)\n",
    "    return f\n",
    "\n",
    "load_csv = os.path.join(path, \"Loads.csv\")\n",
    "\n",
    "with open (load_csv,'r') as f:\n",
    "    f = remove_comments(f)\n",
    "    loads = pd.read_csv(f)\n",
    "    f.close()\n",
    "\n",
    "load_shapes = os.path.join(path, \"LoadShapes.csv\")    \n",
    "    \n",
    "with open (load_shapes,'r') as f:\n",
    "    f = remove_comments(f)\n",
    "    loadshapes = pd.read_csv(f)\n",
    "    f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandapower.plotting as plot\n",
    "net = ieee_european_lv_asymmetric('on_peak_566')\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "try:\n",
    "    import seaborn\n",
    "    colors = seaborn.color_palette()\n",
    "except:\n",
    "    colors = [\"b\", \"g\", \"r\", \"c\", \"y\"]\n",
    "\n",
    "sizes = plot.get_collection_sizes(net)\n",
    "\n",
    "bc = plot.create_bus_collection(net, net.bus.index, size=sizes['bus'], color=colors[0], zorder=10)\n",
    "tlc, tpc = plot.create_trafo_collection(net, net.trafo.index, color=\"c\",size=sizes['trafo'])\n",
    "lcd = plot.create_line_collection(net, net.line.index, color=\"grey\", linewidths=0.5, use_bus_geodata=True)\n",
    "sc = plot.create_ext_grid_collection(net, ext_grid_buses=net.ext_grid.bus.values, size=sizes['ext_grid'], color=\"c\", zorder=11)\n",
    "ldA = plot.create_bus_collection(net, net.asymmetric_load.bus.values[np.where(net.asymmetric_load.p_a_mw >0)], patch_type=\"poly3\", size=sizes['bus'], color=\"r\", zorder=11)\n",
    "ldB = plot.create_bus_collection(net, net.asymmetric_load.bus.values[np.where(net.asymmetric_load.p_b_mw >0)], patch_type=\"rect\", size=sizes['bus'], color=\"y\", zorder=11)\n",
    "ldC = plot.create_bus_collection(net, net.asymmetric_load.bus.values[np.where(net.asymmetric_load.p_c_mw >0)], patch_type=\"circle\", size=sizes['bus'], color=\"g\", zorder=11)\n",
    "plot.draw_collections([lcd, bc, tlc, tpc, sc,ldA,ldB,ldC], figsize=(10,7))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "- Cyan rectangle is External Grid and circle is the transformer\n",
    "- Red triangles, yellow rectangles and green dots represent A,B,C phase loads respectively\n",
    "- blue dots are buses\n",
    "- grey lines are power lines.\n",
    "\n",
    "\n",
    "# Time-series calculation\n",
    "\n",
    "Function to load time series data and perform time-series **three phase load flow**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def timeseries_example(output_dir):\n",
    "    # 1. create test net\n",
    "    net = ieee_european_lv_asymmetric('on_peak_566')\n",
    "    runpp_3ph(net)\n",
    "    # 2. create data source for loads\n",
    "    profiles, ds = create_data_source()\n",
    "    # 3. create controllers (to control P values of the load and the sgen)\n",
    "    net = create_controllers(net, ds)\n",
    "\n",
    "    # time steps to be calculated. Could also be a list with non-consecutive time steps\n",
    "    time_steps = range(0, 672)\n",
    "    print(time_steps)\n",
    "\n",
    "    # 4. the output writer with the desired results to be stored to files.\n",
    "    ow = create_output_writer(net, time_steps, output_dir)\n",
    "\n",
    "    # 5. the main time series function\n",
    "    run_timeseries(net, time_steps, output_writer=ow, run=runpp_3ph, continue_on_divergence=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data source\n",
    "Data taken from csv data provided in IEEE benchmark grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_data_source():\n",
    "    profiles = pd.DataFrame()\n",
    "    for loadprofile,file in (loadshapes[['Name','File']].values):\n",
    "        file_path = os.path.join(load_path, file)\n",
    "        profiles[loadprofile] = pd.read_csv(file_path).mult.values * 1e-3 \n",
    "        profiles[loadprofile+'_Q'] = profiles[loadprofile] * np.tan(\n",
    "                    np.arccos(float(loads[loads.Yearly==loadprofile].PF.values)))    \n",
    "    ds = DFData(profiles)\n",
    "    return profiles, ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Time series controller\n",
    "P, Q values entered using P and power factor(cos_phi) data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_controllers(net, ds):\n",
    "    ConstControl(net, element='asymmetric_load', variable='p_a_mw', element_index=loads[loads['phases']=='A'].index,\n",
    "                 data_source=ds, profile_name=loads[loads['phases']=='A'].Yearly)\n",
    "    ConstControl(net, element='asymmetric_load', variable='q_a_mvar', element_index=loads[loads['phases']=='A'].index,\n",
    "                 data_source=ds, profile_name=loads[loads['phases']=='A'].Yearly+'_Q')\n",
    "    \n",
    "    ConstControl(net, element='asymmetric_load', variable='p_b_mw', element_index=loads[loads['phases']=='B'].index,\n",
    "                 data_source=ds, profile_name=loads[loads['phases']=='B'].Yearly)\n",
    "    ConstControl(net, element='asymmetric_load', variable='q_b_mvar', element_index=loads[loads['phases']=='B'].index,\n",
    "                 data_source=ds, profile_name=loads[loads['phases']=='B'].Yearly+'_Q')\n",
    "    \n",
    "    ConstControl(net, element='asymmetric_load', variable='p_c_mw', element_index=loads[loads['phases']=='C'].index,\n",
    "                 data_source=ds, profile_name=loads[loads['phases']=='C'].Yearly)\n",
    "    ConstControl(net, element='asymmetric_load', variable='q_c_mvar', element_index=loads[loads['phases']=='C'].index,\n",
    "                 data_source=ds, profile_name=loads[loads['phases']=='C'].Yearly+'_Q')   \n",
    "    return net\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Format output\n",
    "- Maximum bus voltage magnitude for each time step\n",
    "- Aggregated Sum of P & Q values for each time step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_output_writer(net, time_steps, output_dir):\n",
    "    ow = OutputWriter(net, time_steps, output_path=output_dir, output_file_type=\".json\")\n",
    "    ow.log_variable('res_trafo_3ph', 'p_a_lv_mw', index=net.trafo.index)\n",
    "    ow.log_variable('res_trafo_3ph', 'p_b_lv_mw', index=net.trafo.index)\n",
    "    ow.log_variable('res_trafo_3ph', 'p_c_lv_mw', index=net.trafo.index)\n",
    "    ow.log_variable('res_trafo_3ph', 'q_a_lv_mvar', index=net.trafo.index)\n",
    "    ow.log_variable('res_trafo_3ph', 'q_b_lv_mvar', index=net.trafo.index)\n",
    "    ow.log_variable('res_trafo_3ph', 'q_c_lv_mvar', index=net.trafo.index)\n",
    "    ow.log_variable('res_bus_3ph', 'vm_a_pu', index=[34])   # Load 1 Phase A \n",
    "    ow.log_variable('res_bus_3ph', 'vm_b_pu', index=[899])  # Load 53 Phase B\n",
    "    ow.log_variable('res_bus_3ph', 'vm_c_pu', index=[614])  # Load 32 Phase C\n",
    "    ow.log_variable('res_bus_3ph', 'va_a_degree', index=[34])   # Load 1 Phase A \n",
    "    ow.log_variable('res_bus_3ph', 'va_b_degree', index=[899])  # Load 53 Phase B\n",
    "    ow.log_variable('res_bus_3ph', 'va_c_degree', index=[614])  # Load 32 Phase C\n",
    "    return ow"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Main Function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_dir = os.path.join(tempfile.gettempdir(), \"time_series_example\")\n",
    "print(\"Results can be found in your local temp folder: {}\".format(output_dir))\n",
    "if not os.path.exists(output_dir):\n",
    "    os.mkdir(output_dir)\n",
    "timeseries_example(output_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Importing data for plotting\n",
    "A pandas dataframe is made from the output '.json' file\n",
    "This dataframe will be used for plotting results using matplotlib\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_dir = os.path.join(tempfile.gettempdir(), \"time_series_example\")\n",
    "# Power Values from secondary of the transformer\n",
    "PA = os.path.join(output_dir,'res_trafo_3ph','p_a_lv_mw.json')\n",
    "PB = os.path.join(output_dir,'res_trafo_3ph','p_b_lv_mw.json')\n",
    "PC = os.path.join(output_dir,'res_trafo_3ph','p_c_lv_mw.json')\n",
    "QA = os.path.join(output_dir,'res_trafo_3ph','q_a_lv_mvar.json')\n",
    "QB = os.path.join(output_dir,'res_trafo_3ph','q_b_lv_mvar.json')\n",
    "QC = os.path.join(output_dir,'res_trafo_3ph','q_c_lv_mvar.json')\n",
    "\n",
    "#pandapower Results\n",
    "df_pp = pd.read_json(PA)*-1e3\n",
    "df_pp['PB'] = pd.read_json(PB)*-1e3\n",
    "df_pp['PC'] = pd.read_json(PC)*-1e3\n",
    "df_pp['P_SUM'] =  df_pp.sum(axis=1)\n",
    "df_pp = df_pp.rename(columns={0:'PA'}) \n",
    "\n",
    "df_pq = pd.read_json(QA)*-1e3\n",
    "df_pq['QB'] = pd.read_json(QB)*-1e3\n",
    "df_pq['QC'] = pd.read_json(QC)*-1e3\n",
    "df_pq = df_pq.rename(columns={0:'QA'})\n",
    "df_pp['Q_SUM'] =  df_pq.sum(axis=1)\n",
    "\n",
    "\n",
    "#The magnitude of voltage at LOAD1 (phase A), LOAD32 (phase C), and LOAD53 (phase B) over the one-day period are shown \n",
    "pp_va = os.path.join(output_dir,'res_bus_3ph','vm_a_pu.json')\n",
    "pp_vb = os.path.join(output_dir,'res_bus_3ph','vm_b_pu.json')\n",
    "pp_vc = os.path.join(output_dir,'res_bus_3ph','vm_c_pu.json')\n",
    "\n",
    "#pandapower Results\n",
    "df_pp_v = pd.read_json(pp_va)*(416/np.sqrt(3))\n",
    "df_pp_v['VB'] = pd.read_json(pp_vb)*(416/np.sqrt(3))\n",
    "df_pp_v['VC'] = pd.read_json(pp_vc)*(416/np.sqrt(3))\n",
    "df_pp_v= df_pp_v.rename(columns={34: \"VA\"})\n",
    "\n",
    "# This is required since json makes keys as string type, the index order is like 1, 10, 100 ,...\n",
    "df_pp_v.index = df_pp_v.index.astype(np.int64)\n",
    "df_pp_v = df_pp_v.sort_index()\n",
    "\n",
    "df_pp.index = df_pp.index.astype(np.int64)\n",
    "df_pp = df_pp.sort_index()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plotting Results\n",
    "\n",
    "The following data is generated :\n",
    "- pandapower voltage and power (time-series values)\n",
    "- OpenDSS voltage and power (not executed here)\n",
    "- Difference between pandapower OpenDSS (time-series values)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "fig,ax1 = plt.subplots(1,1)\n",
    "df_pp_v.plot(kind='line', y=['VA','VB','VC'], ax=ax1, color=['blue','orange','green'],figsize=(10,5),legend = True)\n",
    "ax1.set_ylim((235,255))\n",
    "ax1.xaxis.set_label_text(\"time step\")\n",
    "ax1.yaxis.set_label_text(\"voltage mag. [Volts]\")\n",
    "ax1.set_title(\"Load Voltages (pandapower) / V\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Comparison with OpenDSS Solutions\n",
    "\n",
    "**Voltage values are taken from three loads of each phase :**\n",
    "\n",
    "Load #1 Phase A\n",
    "\n",
    "Load #53 Phase B\n",
    "\n",
    "Load #32 Phase C\n",
    "\n",
    "\n",
    "\n",
    "<img src=\"pics/OpenDSS_time_series.png\" width=\"50%\">\n",
    "\n",
    "<img src=\"pics/dss_fig_Vm.png\" width=\"55%\">"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
